O objetivo desta atividade é aplicar as técnicas de regressão linear. A atividade é dividida em três partes:
Algumas recomendações:
foo não está instalado. Para instalar um pacote, execute o comando install.packages("foo") no Console, ou clique em Tools -> Install Packages.Nesta atividade, nenhuma configuração é necessária por padrão, mas você pode usar o bloco abaixo para incluir o que julgar necessário.
# insira seus comandos de configuração aqui
Nesta parte você irá exercitar a sua familiaridade com o coeficiente de correlação de Pearson, que é uma medida da correlação linear entre duas variáveis. Para isso, basta jogar Guess the Correlation, um jogo online em que você faz pontos ao adivinhar de forma aproximada o coeficiente de correlação para um dado gráfico de dispersão. O jogo está disponível em http://guessthecorrelation.com/. Tire uma selfie onde apareçam seu rosto e seu resultado no jogo (só vale escore de 25 ou mais!), e insira-a no lugar do meu resultado de 59 pontos abaixo:
Guess the Correlation
38pts!
Uma etapa preliminar em uma análise de regressão é estabelecer se as variáveis consideradas efetivamente possuem a relação esperada pelo modelo de regressão. Para isso, utiliza-se o ferramental da análise exploratória de dados.
O arquivo variaveis.dat possui cinco pares de variáveis (x1-y1, x2-y2, …, x5-y5). Nesta parte da atividade, você deve determinar para quais desses pares parece possível construir um modelo de regressão linear. Use o bloco R abaixo para realizar as análises de dados, e nas respostas indique, para cada par x-y, se seria viável uma regressão ou não, justificando sua resposta. Algumas observações:
pares <- read.table("variaveis.dat", header = TRUE)
Respostas aqui
plot(pares$x1, pares$y1)
(cor.x1y1 = cor(pares$x1, pares$y1))
## [1] 0.01115994
(test.x1y1 = cor.test(pares$x1, pares$y1))
##
## Pearson's product-moment correlation
##
## data: pares$x1 and pares$y1
## t = 0.077323, df = 48, p-value = 0.9387
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2680203 0.2886111
## sample estimates:
## cor
## 0.01115994
#summary(test.x1y1)
plot(pares$x2, pares$y2)
(cor.x2y2 = cor(pares$x2, pares$y2))
## [1] 0.9857316
(test.x2y2 = cor.test(pares$x2, pares$y2))
##
## Pearson's product-moment correlation
##
## data: pares$x2 and pares$y2
## t = 40.573, df = 48, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.9748631 0.9919201
## sample estimates:
## cor
## 0.9857316
plot(pares$x3, pares$y3)
(cor.x3y3 = cor(pares$x3, pares$y3))
## [1] 0.7428404
(test.x3y3 = cor.test(pares$x3, pares$y3))
##
## Pearson's product-moment correlation
##
## data: pares$x3 and pares$y3
## t = 7.6875, df = 48, p-value = 6.553e-10
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.5855698 0.8462175
## sample estimates:
## cor
## 0.7428404
#sd(pares$x3)
#sd(pares$y3)
plot(pares$x4, pares$y4)
(cor.x4y4 = cor(pares$x4, pares$y4))
## [1] -0.972956
(test.x4y4 = cor.test(pares$x4, pares$y4))
##
## Pearson's product-moment correlation
##
## data: pares$x4 and pares$y4
## t = -29.182, df = 48, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.9846427 -0.9525884
## sample estimates:
## cor
## -0.972956
plot(pares$x5, pares$y5)
(cor.x5y5 = cor(pares$x5, pares$y5))
## [1] 0.07415827
(test.x5y5 = cor.test(pares$x5, pares$y5))
##
## Pearson's product-moment correlation
##
## data: pares$x5 and pares$y5
## t = 0.5152, df = 48, p-value = 0.6088
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.2084931 0.3453768
## sample estimates:
## cor
## 0.07415827
Para essa comparação iremos usar tempos de execução medidos pelo script Python fibo.py. Esse script mede o tempo necessário para encontrar N vezes o quinquagésimo termo da série de Fibonacci. O script varia o valor de N de acordo com três parâmetros na linha de comando: ~ - o valor inicial de N; - o valor final de N; - o passo de incremento.
O quarto parâmetro na linha de comando é quantas repetições são realizadas para cada valor testado; se o parâmetro não for especificado, cada valor de N é mensurado apenas uma vez. Por exemplo, o comando
python3 fibo.py 100 200 25 2
1200000 - 1000000 - 120000 executa o script com N valendo 100, 125, 150, 175 e 200, com duas repetições para cada valor de N.
O script pode ser executado no RStudio Cloud. Na janela inferior esquerda, normalmente usada para o console, há uma aba Terminal, na qual você pode executar comandos do Linux.
Os passos deste experimento são os seguintes:
Execute o script fibo.py. Os parâmetros devem ser ajustados a seu critério, respeitando as seguintes considerações:
Especifique um modelo de regressão para o tempo de execução em função de N e ajuste-o aos dados coletados no passo 1. Qual a equação de regressão obtida?
Verifique a qualidade do ajuste e faça o diagnóstico dos resíduos do modelo. Determine se o modelo é adequado ou não.
Use o modelo para obter predições para quatro valores de N diferentes dos usados no script, sendo dois valores dentro do intervalo mensurado para N e dois valores fora desse intervalo (por exemplo, considerando o exemplo de execução do script, onde o intervalo era 100-200, os valores poderiam ser 51, 123, 171 e 500). Além das estimativas pontuais, obtenha também os intervalos de confiança de 95% para as predições.
Execute novamente o script para obter uma medida para cada um dos valores de N usados no item anterior. Calcule os erros \(e_i=y_i-\hat{y}_i\) correspondentes a essas medidas, e verifique se as medidas estão dentro dos intervalos de confiança obtidos no item 4.
De acordo com o modelo, qual o valor de N correspondente a um tempo de execução de 23 segundos? Execute o script para obter o tempo de execução para esse valor de N, e determine o erro \(e_i\) dessa medida.
Lembre-se que os dados experimentais devem ser salvos em arquivos para que sua análise seja reproduzível. Para facilitar essa tarefa, o script já gera a saída em um formato apropriado; você pode redirecionar a saída do script para um arquivo (por exemplo, python3 fibo.py 100 200 25 2 >fibo.dat) ou simplesmente criar o arquivo de dados no próprio editor do RStudio Cloud (crie um novo arquivo texto e cole a saída do script). Para as execuções dos passos 5 e 6 você pode gerar um único arquivo de dados.
# seu código R aqui
# python3 fibo.py 120000 1000000 5000 5 > part3.dat
# python3 fibo.py 90000 1300000 5000 5 > part3_12.dat
part3 <- read.table("part3.dat", header = TRUE)
#length(part3)
str(part3)
## 'data.frame': 885 obs. of 2 variables:
## $ n : int 120000 120000 120000 120000 120000 125000 125000 125000 125000 125000 ...
## $ tempo: num 0.58 0.579 0.585 0.584 0.584 ...
lmp3 = lm(tempo ~ n, data=part3)
lmp3
##
## Call:
## lm(formula = tempo ~ n, data = part3)
##
## Coefficients:
## (Intercept) n
## -1.712e-02 4.746e-06
#str(lmp3)
(lmp3.summary = summary(lmp3))
##
## Call:
## lm(formula = tempo ~ n, data = part3)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.141692 -0.033848 -0.007993 0.033530 0.144072
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -1.712e-02 3.875e-03 -4.417 1.12e-05 ***
## n 4.746e-06 6.296e-09 753.811 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.04785 on 883 degrees of freedom
## Multiple R-squared: 0.9984, Adjusted R-squared: 0.9984
## F-statistic: 5.682e+05 on 1 and 883 DF, p-value: < 2.2e-16
#str(lmp3.summary)
#par(mfrow=c(2,2));
plot(lmp3)
plot(part3)
abline(lmp3, col='red')
diffResidxFitted = 100*resid(lmp3)/fitted(lmp3)
max(diffResidxFitted)
## [1] 5.976154
(lmp3.predict = predict(lmp3, newdata = data.frame(n=c(100000, 400000, 800000, 1200000)), conf.level=0.95))
## 1 2 3 4
## 0.4574549 1.8811743 3.7794667 5.6777591
(tempo.400000 = part3$tempo[part3$n == 400000])
## [1] 1.84171 1.84331 1.84262 1.84096 1.84353
#tempo.400000.diff = diff(tempo.400000, lmp3.predict[2])
tempo.400000.diff = abs(tempo.400000 - lmp3.predict[2])
(tempo.800000 = part3$tempo[part3$n == 800000])
## [1] 3.84203 3.84017 3.85909 3.86171 3.85011
#tempo.800000.diff = diff(tempo.800000, lmp3.predict[3])
tempo.800000.diff = abs(tempo.800000 - lmp3.predict[3])
part312 <- read.table("part3_12.dat", header = TRUE)
(tempo.100000 = part312$tempo[part312$n == 100000])
## [1] 0.46315 0.46320 0.46335 0.46445 0.46355
tempo.100000.diff = abs(tempo.100000 - lmp3.predict[1])
(tempo.1200000 = part312$tempo[part312$n == 1200000])
## [1] 5.56335 5.56422 5.56622 5.57196 5.56985
tempo.1200000.diff = abs(tempo.1200000 - lmp3.predict[4])
(lmp3.predict.confinterval = lmp3.predict - (lmp3.predict * 0.95))
## 1 2 3 4
## 0.02287275 0.09405871 0.18897333 0.28388796
lmp3.predict[4]
## 4
## 5.677759
tempo.100000.diff < lmp3.predict.confinterval[1]
## [1] TRUE TRUE TRUE TRUE TRUE
tempo.400000.diff < lmp3.predict.confinterval[2]
## [1] TRUE TRUE TRUE TRUE TRUE
tempo.800000.diff < lmp3.predict.confinterval[3]
## [1] TRUE TRUE TRUE TRUE TRUE
tempo.1200000.diff < lmp3.predict.confinterval[4]
## [1] TRUE TRUE TRUE TRUE TRUE
(tempo.23s = (lmp3.summary$coefficients[1] - 23) / (-lmp3.summary$coefficients[2]))
## [1] 4850068
predict(lmp3, newdata = data.frame(n=c(tempo.23s)), conf.level=0.95)
## 1
## 23
tempo.23s
## [1] 4850068
#tempo.23s =
part23s <- read.table("part3_4850068.dat", header = TRUE)
part23s
## n tempo
## 1 4850068 23.14015
## 2 4850068 23.35855
## 3 4850068 23.34468
## 4 4850068 23.35834
## 5 4850068 23.36985
tempo.23s
## [1] 4850068
(part23s.conf = (23 * 100) / part23s$tempo)
## [1] 99.39434 98.46502 98.52352 98.46590 98.41741
(part23s.err = 100 - part23s.conf)
## [1] 0.6056573 1.5349840 1.4764820 1.5340987 1.5825947
Respostas aqui 1 - CLI: 120000 inicia em 0.54 e 1200000 termina em 5.4 - python3 fibo.py 120000 1000000 5000 5 > part3.dat
2 - Formula = tempo ~ n. Equaçao de regressão obtida: tempo = -0.0171 + 5^{-6} n
3 - O resíduo do modelo, está meio disperso como se fosse um cone pelo apresentado no gráfico, contudo ao mensurar a máxima diferença entre os resíduos e ajustes chega só a 5.98 %, e os resíduos na normal, ficou bacana como se fosse um normal ou no máximo ‘light tailed’ com pequenas mudanças nas caudas, o R-squared é bom e o modelo é significativo então meu veredito final é que o modelo é adequado.
4: N x Estimativa
100000 x 0.457s, Nível de Confiança 95%: 0.023s
400000 x 1.881s, Nível de Confiança 95%: 0.094s
800000 x 3.779s, Nível de Confiança 95%: 0.189s
1200000 x 5.678s, Nível de Confiança 95%: 0.284s
Comparando as duas predições (400000, 800000) com as medicoes executadas
1.881 - 1.842, 1.843, 1.843, 1.841, 1.844 = 0.039, 0.038, 0.039, 0.04, 0.038
3.779 - 3.842, 3.84, 3.859, 3.862, 3.85 = 0.063, 0.061, 0.08, 0.082, 0.071
Devo reconhecer que eu esperava uma diferença menor da predição que ja existia nos dados utilizados para especificação do modelo, mas em todos os dois casos, a predição fica dentro do intervalo de confiança de 95%
5: Acredito que o modelo se adaptou bem para os valores que não estavam contemplados nos dados coletados e todos foram validos para um Nível de confiança de 95%, porem se fosse para 99% nenhuma prediçao seria valida
100000: 0.457s 0.463, 0.463, 0.463, 0.464, 0.464 = 0.006, 0.006, 0.006, 0.007, 0.006 < Nível de Confiança 95%: TRUE, TRUE, TRUE, TRUE, TRUE
1200000: 5.678s - 5.563, 5.564, 5.566, 5.572, 5.57 = 0.114, 0.114, 0.112, 0.106, 0.108 < Nível de Confiança 95%: TRUE, TRUE, TRUE, TRUE, TRUE
6: De acordo com o modelo, o valor de N correspondente a 23s é de 4.8500679^{6}, Ja testando o mesmo N 4.8500679^{6} no script python3 fibo.py 4850068 4850068 5000 5, gera os tempos de 23.14015, 23.35855, 23.34468, 23.35834, 23.36985 s. E os erros relacionado ao 23s são de 0.61, 1.53, 1.48, 1.53, 1.58 % respectivamente. O que é considerado algo bom utilizando um Nível de confiança de 95%, porém invalido para algumas medidas utilizando 99%.